3. Comparison of IsoHD splicing events with post-mortem HD patient data

In [1]:
library(dplyr)
library(ggplot2)
library(biomaRt)
library(org.Hs.eg.db)
library(clusterProfiler)
library(pheatmap)
library(DT)
library(dendextend)
library(dendsort)
library(ggplotify)
library(gridExtra)
library(reshape2)
library(Mfuzz)
library(enrichplot)
Attaching package: ‘dplyr’


The following objects are masked from ‘package:stats’:

    filter, lag


The following objects are masked from ‘package:base’:

    intersect, setdiff, setequal, union


Warning message:
“package ‘ggplot2’ was built under R version 4.0.5”
Loading required package: AnnotationDbi

Loading required package: stats4

Loading required package: BiocGenerics

Loading required package: parallel


Attaching package: ‘BiocGenerics’


The following objects are masked from ‘package:parallel’:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB


The following objects are masked from ‘package:dplyr’:

    combine, intersect, setdiff, union


The following objects are masked from ‘package:stats’:

    IQR, mad, sd, var, xtabs


The following objects are masked from ‘package:base’:

    anyDuplicated, append, as.data.frame, basename, cbind, colnames,
    dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
    grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
    order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
    rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
    union, unique, unsplit, which.max, which.min


Loading required package: Biobase

Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.


Loading required package: IRanges

Loading required package: S4Vectors


Attaching package: ‘S4Vectors’


The following objects are masked from ‘package:dplyr’:

    first, rename


The following object is masked from ‘package:base’:

    expand.grid



Attaching package: ‘IRanges’


The following objects are masked from ‘package:dplyr’:

    collapse, desc, slice



Attaching package: ‘AnnotationDbi’


The following object is masked from ‘package:dplyr’:

    select






clusterProfiler v3.18.0  For help: https://guangchuangyu.github.io/software/clusterProfiler

If you use clusterProfiler in published research, please cite:
Guangchuang Yu, Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.


Attaching package: ‘clusterProfiler’


The following object is masked from ‘package:AnnotationDbi’:

    select


The following object is masked from ‘package:IRanges’:

    slice


The following object is masked from ‘package:S4Vectors’:

    rename


The following object is masked from ‘package:biomaRt’:

    select


The following object is masked from ‘package:stats’:

    filter



---------------------
Welcome to dendextend version 1.15.1
Type citation('dendextend') for how to cite the package.

Type browseVignettes(package = 'dendextend') for the package vignette.
The github page is: https://github.com/talgalili/dendextend/

Suggestions and bug-reports can be submitted at: https://github.com/talgalili/dendextend/issues
Or contact: <tal.galili@gmail.com>

	To suppress this message use:  suppressPackageStartupMessages(library(dendextend))
---------------------



Attaching package: ‘dendextend’


The following object is masked from ‘package:stats’:

    cutree


Warning message:
“package ‘gridExtra’ was built under R version 4.0.5”

Attaching package: ‘gridExtra’


The following object is masked from ‘package:Biobase’:

    combine


The following object is masked from ‘package:BiocGenerics’:

    combine


The following object is masked from ‘package:dplyr’:

    combine


Warning message:
“package ‘reshape2’ was built under R version 4.0.5”
Loading required package: e1071

Warning message:
“no DISPLAY variable so Tk is not available”

Attaching package: ‘widgetTools’


The following object is masked from ‘package:dplyr’:

    funs



Attaching package: ‘DynDoc’


The following object is masked from ‘package:BiocGenerics’:

    path


3.1. VastDB annotated splicing events in IsoHD, HD cortex, HD striatum

In [2]:
leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv <- read.table(
    file='leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv',
    sep='\t'
)

leafcutterAStable_human_cortex_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv <- read.table(
    file='leafcutterAStable_human_cortex_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv',
    sep='\t'
)

leafcutterAStable_human_striatum_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv <- read.table(
    file='leafcutterAStable_human_striatum_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv',
    sep='\t'
)
In [3]:
options(repr.plot.height=6,repr.plot.width=6)
ggvenn::ggvenn(
    data=list(
        'IsoHD'=leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        subset(DiffSplicing.test!='Cell_type_NPCvsESC') %>% 
        with(gsub('-[0-9]/[0-9]$','',EVENT) %>% na.omit %>% unique),
        'Cortex HD'=leafcutterAStable_human_cortex_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        with(gsub('-[0-9]/[0-9]$','',EVENT) %>% na.omit %>% unique),
        'Striatum HD'=leafcutterAStable_human_striatum_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        with(gsub('-[0-9]/[0-9]$','',EVENT) %>% na.omit %>% unique)
    ),
    fill_color=c("blue","yellow","green", "red"),
    show_percentage=FALSE,
    set_name_size=10,
    text_size=10
)
In [4]:
list(
    'IsoHD'=leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    subset(DiffSplicing.test!='Cell_type_NPCvsESC') %>% 
    with(JunctionName %>% na.omit %>% unique %>% sort),
    'Cortex HD'=leafcutterAStable_human_cortex_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    with(JunctionName %>% na.omit %>% unique %>% sort),
    'Striatum HD'=leafcutterAStable_human_striatum_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    with(JunctionName %>% na.omit %>% unique %>% sort)
) %>% 
lapply(length)
$IsoHD
5255
$`Cortex HD`
332
$`Striatum HD`
5919
In [5]:
leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents.tsv %>% 
subset(DiffSplicing.test!='Cell_type_NPCvsESC' & (EVENT %in% (Reduce(intersect,list(
    'IsoHD'=leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    subset(DiffSplicing.test!='Cell_type_NPCvsESC') %>% 
    with(EVENT %>% na.omit %>% unique),
    'Cortex HD'=leafcutterAStable_human_cortex_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    with(EVENT %>% na.omit %>% unique),
    'Striatum HD'=leafcutterAStable_human_striatum_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    with(EVENT %>% na.omit %>% unique)
))))) %>% 
with(genes %>% unique %>% sort)
  1. 'BCOR'
  2. 'CAMK2G'
  3. 'CPSF7'
  4. 'EHBP1'
  5. 'KCNMA1'
  6. 'PTPRM'
  7. 'RPRD2'
  8. 'TBC1D5'

3.2. Overlapping IsoHD-cortex HD-striatum HD DSJs

In [6]:
intersect(
    leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    subset(grepl('_disease_stage_',DiffSplicing.test)) %>% 
    with(gsub('-[0-9]/[0-9]$','',EVENT)) %>% na.omit %>% unique,
    leafcutterAStable_human_striatum_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    with(gsub('-[0-9]/[0-9]$','',EVENT)) %>% na.omit %>% unique
) %>% unique %>% length

intersect(
    leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    subset(grepl('_disease_stage_',DiffSplicing.test)) %>% 
    with(gsub('-[0-9]/[0-9]$','',EVENT)) %>% na.omit %>% unique,
    leafcutterAStable_human_cortex_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
    with(gsub('-[0-9]/[0-9]$','',EVENT)) %>% na.omit %>% unique
) %>% unique %>% length
280
14
In [7]:
intersect(
        leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        subset(grepl('_disease_stage_',DiffSplicing.test)) %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique,
        leafcutterAStable_human_striatum_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique
    ) %>% unique %>% length

intersect(
        leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        subset(grepl('_disease_stage_',DiffSplicing.test)) %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique,
        leafcutterAStable_human_cortex_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique
    ) %>% unique %>% length
988
53
In [8]:
c(
    intersect(
        leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        subset(grepl('_disease_stage_',DiffSplicing.test)) %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique,
        leafcutterAStable_human_striatum_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique
    ),
    intersect(
        leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        subset(grepl('_disease_stage_',DiffSplicing.test)) %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique,
        leafcutterAStable_human_cortex_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique
    )
) %>% unique %>% length
1012
In [9]:
isohdcortexstriatum.alljunctions.naomit <- c(
    intersect(
        leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        subset(grepl('_disease_stage_',DiffSplicing.test)) %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique,
        leafcutterAStable_human_striatum_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique
    ),
    intersect(
        leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        subset(grepl('_disease_stage_',DiffSplicing.test)) %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique,
        leafcutterAStable_human_cortex_HD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
        with(gsub('clu_[0-9]*_','',JunctionName)) %>% unique
    )
) %>% unique %>% 
(function(juncnames) {
    Reduce(f=intersect,x=list(
        readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
            juncnames,
            coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('hESC',Celltype)) %>% row.names
        ],
        readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
            juncnames,
            coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('NPC',Celltype)) %>% row.names
        ],
        readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
            juncnames,
            coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('Neuron',Celltype)) %>% row.names
        ],
        readsunique.novel.human_cortex_hd_perind_logit.constcounts.known2[juncnames,],
        readsunique.novel.human_striatum_hd_perind_logit.constcounts.known2[juncnames,]
    ) %>% 
    lapply(function(x) {
        x %>% 
        t %>% scale %>% t %>% 
        na.omit %>% 
        row.names
    }))
}) %>% 
unlist %>% unique
In [10]:
isohdcortexstriatum.alljunctions.naomit %>% length
710

3.3. Fuzzy clustering of shared DSJs

In [11]:
isohdcortexstriatum.Dmin <- cbind(
    readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
        isohdcortexstriatum.alljunctions.naomit,
        coldata_ESCNPCNeu_isoHD2 %>% subset(Celltype=='hESC') %>% row.names
    ] %>% 
    t %>% scale %>% t,
    readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
        isohdcortexstriatum.alljunctions.naomit,
        coldata_ESCNPCNeu_isoHD2 %>% subset(Celltype=='NPC') %>% row.names
    ] %>% 
    t %>% scale %>% t,
    readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
        isohdcortexstriatum.alljunctions.naomit,
        coldata_ESCNPCNeu_isoHD2 %>% subset(Celltype=='Neuron') %>% row.names
    ] %>% 
    t %>% scale %>% t
) %>% 
ExpressionSet(
    phenoData=coldata_ESCNPCNeu_isoHD2[
        c(coldata_ESCNPCNeu_isoHD2 %>% subset(Celltype=='hESC') %>% row.names,coldata_ESCNPCNeu_isoHD2[-c(8,14,21),] %>% subset(Celltype=='NPC') %>% row.names,coldata_ESCNPCNeu_isoHD2[-c(8,14,21),] %>% subset(Celltype=='Neuron') %>% row.names),
    ] %>% AnnotatedDataFrame
) %>% 
(function(myexprset) {
    myexprset %>% 
    Dmin(
        m=myexprset %>% mestimate,
        crange=2:10,
        repeats=5,
        visu=FALSE
    )
})
In [12]:
options(repr.plot.width=10,repr.plot.height=5)
ggplot(
    mapping=aes(
        x=2:10,
        y=isohdcortexstriatum.Dmin
    )
) +
geom_point() +
geom_line() +
geom_vline(xintercept=4,color='red',linetype='dashed') +
labs(x='Num of clusters',y='Dmin') +
theme_bw(base_size=16)
In [13]:
isohdcortexstriatum.mfuzz.list <- c(2:5) %>% 
lapply(function(x) {
    myexprset <- cbind(
        readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
            isohdcortexstriatum.alljunctions.naomit,
            coldata_ESCNPCNeu_isoHD2 %>% subset(Celltype=='hESC') %>% row.names
        ] %>% 
        t %>% scale %>% t,
        readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
            isohdcortexstriatum.alljunctions.naomit,
            coldata_ESCNPCNeu_isoHD2 %>% subset(Celltype=='NPC') %>% row.names
        ] %>% 
        t %>% scale %>% t,
        readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
            isohdcortexstriatum.alljunctions.naomit,
            coldata_ESCNPCNeu_isoHD2 %>% subset(Celltype=='Neuron') %>% row.names
        ] %>% 
        t %>% scale %>% t
    ) %>% 
    ExpressionSet(
        phenoData=coldata_ESCNPCNeu_isoHD2[
            c(coldata_ESCNPCNeu_isoHD2 %>% subset(Celltype=='hESC') %>% row.names,coldata_ESCNPCNeu_isoHD2[-c(8,14,21),] %>% subset(Celltype=='NPC') %>% row.names,coldata_ESCNPCNeu_isoHD2[-c(8,14,21),] %>% subset(Celltype=='Neuron') %>% row.names),
        ] %>% AnnotatedDataFrame
    )
    
    myexprset %>% 
    mfuzz(centers=x,m=myexprset %>% mestimate)
}) %>% setNames(nm=c(2:5))
In [14]:
boxplotglist <- list()

boxplotglist[['hesc']] <- readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('ESC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
subset(membership>=.6) %>% 
ggplot(
    mapping=aes(
        x=Var2 %>% factor(
            levels=coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('ESC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names,
            labels=gsub('readsunique.novel.|.bed','',coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('ESC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names)
        ),
        y=value
    )
) +
geom_boxplot(outlier.alpha=0) +
geom_hline(yintercept=0,linetype='dashed') +
facet_grid(
    facets=cluster~'isoHD: hESC',
    labeller=labeller(cluster=c('1'='Cluster 1','2'='Cluster 2','3'='Cluster 3','4'='Cluster 4','5'='Cluster 5'))
) +
theme_bw(base_size=16) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5)) +
labs(x='',y='Scaled PSI')

boxplotglist[['npc']] <- readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('NPC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
subset(membership>=.6) %>% 
ggplot(
    mapping=aes(
        x=Var2 %>% factor(
            levels=coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('NPC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names,
            labels=gsub('readsunique.novel.|.bed','',coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('NPC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names)
        ),
        y=value
    )
) +
geom_boxplot(outlier.alpha=0) +
geom_hline(yintercept=0,linetype='dashed') +
facet_grid(
    facets=cluster~'isoHD: NPC',
    labeller=labeller(cluster=c('1'='Cluster 1','2'='Cluster 2','3'='Cluster 3','4'='Cluster 4','5'='Cluster 5'))
) +
theme_bw(base_size=16) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5)) +
labs(x='',y='Scaled PSI')

boxplotglist[['neuron']] <- readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('Neuron',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
subset(membership>=.6) %>% 
ggplot(
    mapping=aes(
        x=Var2 %>% factor(
            levels=coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('Neuron',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names,
            labels=gsub('readsunique.novel.|.bed','',coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('Neuron',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names)
        ),
        y=value
    )
) +
geom_boxplot(outlier.alpha=0) +
geom_hline(yintercept=0,linetype='dashed') +
facet_grid(
    facets=cluster~'isoHD: Neuron',
    labeller=labeller(cluster=c('1'='Cluster 1','2'='Cluster 2','3'='Cluster 3','4'='Cluster 4','5'='Cluster 5'))
) +
theme_bw(base_size=16) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5)) +
labs(x='',y='Scaled PSI')

boxplotglist[['cortex']] <- readsunique.novel.human_cortex_hd_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_cortex_hd %>% arrange(Diagnosis) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
subset(membership>=.6) %>% 
ggplot(
    mapping=aes(
        x=Var2 %>% factor(
            levels=coldata_cortex_hd %>% arrange(Diagnosis) %>% row.names,
            labels=gsub('readsunique.novel.|.bed','',coldata_cortex_hd %>% arrange(Diagnosis) %>% row.names)
        ),
        y=value
    )
) +
geom_boxplot(outlier.alpha=0) +
geom_hline(yintercept=0,linetype='dashed') +
facet_grid(
    facets=cluster~'Cortex HD',
    labeller=labeller(cluster=c('1'='Cluster 1','2'='Cluster 2','3'='Cluster 3','4'='Cluster 4','5'='Cluster 5'))
) +
theme_bw(base_size=16) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5)) +
labs(x='',y='Scaled PSI')

boxplotglist[['striatum']] <- readsunique.novel.human_striatum_hd_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_Striatum_HD %>% arrange(Condition) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
subset(membership>=.6) %>% 
ggplot(
    mapping=aes(
        x=Var2 %>% factor(
            levels=coldata_Striatum_HD %>% arrange(Condition) %>% row.names,
            labels=gsub('readsunique.novel.|.bed','',coldata_Striatum_HD %>% arrange(Condition) %>% row.names)
        ),
        y=value
    )
) +
geom_boxplot(outlier.alpha=0) +
geom_hline(yintercept=0,linetype='dashed') +
facet_grid(
    facets=cluster~'Striatum HD',
    labeller=labeller(cluster=c('1'='Cluster 1: 107','2'='Cluster 2: 129','3'='Cluster 3: 129','4'='Cluster 4: 146','5'='Cluster 5'))
) +
theme_bw(base_size=16) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5)) +
labs(x='',y='Scaled PSI')
In [15]:
options(repr.plot.height=10,repr.plot.width=16)
grid.arrange(grobs=boxplotglist[c(1,2,3,4,5)],layout_matrix=matrix(1:5,nrow=1))
In [16]:
fuzzyglist <- list()

fuzzyglist[['hesc']] <- readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('ESC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
arrange(membership) %>% 
(function(mydf) {
    mydf$Var1 <- mydf$Var1 %>% as.character %>% factor(levels=mydf$Var1 %>% as.character %>% unique)
    mydf
}) %>% 
subset(membership>=.6) %>% 
(function(mydf) {
mydf %>% ggplot(
    mapping=aes(
        x=Var2 %>% factor(
            levels=coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('ESC',Celltype)) %>% arrange(Qlength) %>% row.names,
            labels=gsub('readsunique.novel.|.bed','',coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('ESC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names)
        ),
        y=value,
        group=Var1,
        color=membership
    )
) +
geom_line(alpha=.8) +
geom_line(
    data=mydf %>% 
    group_by(Var2,cluster) %>% 
    summarise(meanvalue=value %>% mean),
    mapping=aes(y=meanvalue,group='mean'),
    color='black',size=2
) +
geom_hline(yintercept=0,linetype='dashed',size=.5) +
scale_color_gradientn(guide=FALSE,
    colours=colorRampPalette(colors=c('yellow','green','blue','purple','red'))(255),
    limits=c(.6,1),
    na.value=NA
) +
facet_grid(
    facets=cluster~'isoHD: hESC',
    labeller=labeller(cluster=c('1'='Cluster 1','2'='Cluster 2','3'='Cluster 3','4'='Cluster 4','5'='Cluster 5'))
) +
theme_bw(base_size=16) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5)) +
labs(x='',y='Scaled PSI')
})

fuzzyglist[['npc']] <- readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('NPC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
arrange(membership) %>% 
(function(mydf) {
    mydf$Var1 <- mydf$Var1 %>% as.character %>% factor(levels=mydf$Var1 %>% as.character %>% unique)
    mydf
}) %>% 
subset(membership>=.6) %>% 
(function(mydf) {
mydf %>% ggplot(
    mapping=aes(
        x=Var2 %>% factor(
            levels=coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('NPC',Celltype)) %>% arrange(Qlength) %>% row.names,
            labels=gsub('readsunique.novel.|.bed','',coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('NPC',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names)
        ),
        y=value,
        group=Var1,
        color=membership
    )
) +
geom_line(alpha=.8) +
geom_line(
    data=mydf %>% 
    group_by(Var2,cluster) %>% 
    summarise(meanvalue=value %>% mean),
    mapping=aes(y=meanvalue,group='mean'),
    color='black',size=2
) +
geom_hline(yintercept=0,linetype='dashed',size=.5) +
scale_color_gradientn(guide=FALSE,
    colours=colorRampPalette(colors=c('yellow','green','blue','purple','red'))(255),
    limits=c(.6,1),
    na.value=NA
) +
facet_grid(
    facets=cluster~'isoHD: NPC',
    labeller=labeller(cluster=c('1'='Cluster 1','2'='Cluster 2','3'='Cluster 3','4'='Cluster 4','5'='Cluster 5'))
) +
theme_bw(base_size=16) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5)) +
labs(x='',y='Scaled PSI')
})

fuzzyglist[['neuron']] <- readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('Neuron',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
arrange(membership) %>% 
(function(mydf) {
    mydf$Var1 <- mydf$Var1 %>% as.character %>% factor(levels=mydf$Var1 %>% as.character %>% unique)
    mydf
}) %>% 
subset(membership>=.6) %>% 
(function(mydf) {
    mydf %>% 
    ggplot(
        mapping=aes(
            x=Var2 %>% factor(
                levels=coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('Neuron',Celltype)) %>% arrange(Qlength) %>% row.names,
                labels=gsub('readsunique.novel.|.bed','',coldata_ESCNPCNeu_isoHD2 %>% subset(grepl('Neuron',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names)
            ),
            y=value,
            group=Var1,
            color=membership
        )
    ) +
    geom_line(alpha=.8) +
    geom_line(
        data=mydf %>% 
        group_by(Var2,cluster) %>% 
        summarise(meanvalue=value %>% mean),
        mapping=aes(y=meanvalue,group='mean'),
        color='black',size=2
    ) +
geom_hline(yintercept=0,linetype='dashed',size=.5) +
    scale_color_gradientn(guide=FALSE,
        colours=colorRampPalette(colors=c('yellow','green','blue','purple','red'))(255),
        limits=c(.6,1),
        na.value=NA
    ) +
    facet_grid(
        facets=cluster~'isoHD: Neuron',
        labeller=labeller(cluster=c('1'='Cluster 1: 107','2'='Cluster 2: 129','3'='Cluster 3: 129','4'='Cluster 4: 146','5'='Cluster 5'))
    ) +
    theme_bw(base_size=16) +
    theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5),legend.title=element_text(size=12),legend.text=element_text(size=10)) +
    labs(x='',y='Scaled PSI')
})

fuzzyglist[['cortex']] <- readsunique.novel.human_cortex_hd_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_cortex_hd %>% arrange(Diagnosis) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
arrange(membership) %>% 
(function(mydf) {
    mydf$Var1 <- mydf$Var1 %>% as.character %>% factor(levels=mydf$Var1 %>% as.character %>% unique)
    mydf
}) %>% 
subset(membership>=.6) %>% 
(function(mydf) {
mydf %>% ggplot(
    mapping=aes(
        x=Var2 %>% factor(
            levels=coldata_cortex_hd %>% arrange(Diagnosis) %>% row.names,
            labels=gsub('readsunique.novel.|.bed','',coldata_cortex_hd %>% arrange(Diagnosis) %>% row.names)
        ),
        y=value,
        group=Var1,
        color=membership
    )
) +
geom_line(color='gray20',alpha=.9) +
geom_line(
    data=mydf %>% 
    group_by(Var2,cluster) %>% 
    summarise(meanvalue=value %>% mean),
    mapping=aes(y=meanvalue,group='mean'),
    color='black',size=2
) +
geom_hline(yintercept=0,linetype='dashed',size=.5) +
scale_color_gradientn(guide=FALSE,
    colours=colorRampPalette(colors=c('yellow','green','blue','purple','red'))(255),
    limits=c(.6,1),
    na.value=NA
) +
facet_grid(
    facets=cluster~'Cortex HD',
    labeller=labeller(cluster=c('1'='Cluster 1','2'='Cluster 2','3'='Cluster 3','4'='Cluster 4','5'='Cluster 5'))
) +
theme_bw(base_size=16) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5)) +
labs(x='',y='Scaled PSI')
})

fuzzyglist[['striatum']] <- readsunique.novel.human_striatum_hd_perind_logit.constcounts.known2[
    isohdcortexstriatum.alljunctions.naomit,
    coldata_Striatum_HD %>% arrange(Condition) %>% row.names
] %>% 
t %>% scale %>% t %>% 
melt %>% 
(function(mydf) {
    mydf %>% 
    data.frame(
        cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
        membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
    )
}) %>% 
arrange(membership) %>% 
(function(mydf) {
    mydf$Var1 <- mydf$Var1 %>% as.character %>% factor(levels=mydf$Var1 %>% as.character %>% unique)
    mydf
}) %>% 
subset(membership>=.6) %>% 
(function(mydf) {
mydf %>% ggplot(
    mapping=aes(
        x=Var2 %>% factor(
            levels=coldata_Striatum_HD %>% arrange(Condition) %>% row.names,
            labels=gsub('readsunique.novel.|.bed','',coldata_Striatum_HD %>% arrange(Condition) %>% row.names)
        ),
        y=value,
        group=Var1,
        color=membership
    )
) +
geom_line(color='gray20',alpha=.8) +
geom_line(
    data=mydf %>% 
    group_by(Var2,cluster) %>% 
    summarise(meanvalue=value %>% mean),
    mapping=aes(y=meanvalue,group='mean'),
    color='black',size=2
) +
geom_hline(yintercept=0,linetype='dashed',size=.5) +
scale_color_gradientn(guide=FALSE,
    colours=colorRampPalette(colors=c('yellow','green','blue','purple','red'))(255),
    limits=c(.6,1),
    na.value=NA
) +
facet_grid(
    facets=cluster~'Striatum HD',
    labeller=labeller(cluster=c('1'='Cluster 1','2'='Cluster 2','3'='Cluster 3','4'='Cluster 4','5'='Cluster 5')),
    drop=FALSE
) +
theme_bw(base_size=16) +
theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5)) +
labs(x='',y='Scaled PSI')
})
`summarise()` regrouping output by 'Var2' (override with `.groups` argument)

`summarise()` regrouping output by 'Var2' (override with `.groups` argument)

`summarise()` regrouping output by 'Var2' (override with `.groups` argument)

`summarise()` regrouping output by 'Var2' (override with `.groups` argument)

`summarise()` regrouping output by 'Var2' (override with `.groups` argument)

In [17]:
options(repr.plot.height=10,repr.plot.width=17)
grid.arrange(grobs=fuzzyglist[c(1,2,3,4,5)],layout_matrix=matrix(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5),nrow=1))
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
In [18]:
options(repr.plot.height=10,repr.plot.width=17)
grid.arrange(grobs=append(fuzzyglist[c(1,2,3)],boxplotglist[c(4,5)]),layout_matrix=matrix(c(1,1,1,2,2,2,3,3,3,4,4,4,5,5,5),nrow=1))
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
In [19]:
options(repr.plot.height=10,repr.plot.width=5)
fuzzylegend <- (
    readsunique.novel.human_ESCNPCNeu_isoHD_perind_logit.constcounts.known2[
        isohdcortexstriatum.alljunctions.naomit,
        coldata_ESCNPCNeu_isoHD2[-c(8,14,21),] %>% subset(grepl('Neuron',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names
    ] %>% 
    t %>% scale %>% t %>% 
    melt %>% 
    (function(mydf) {
        mydf %>% 
        data.frame(
            cluster=isohdcortexstriatum.mfuzz.list[['4']]$cluster[mydf$Var1] %>% factor,
            membership=isohdcortexstriatum.alljunctions.naomit.notstable.membership[mydf$Var1]
        )
    }) %>% 
    arrange(membership) %>% 
    (function(mydf) {
        mydf$Var1 <- mydf$Var1 %>% as.character %>% factor(levels=mydf$Var1 %>% as.character %>% unique)
        mydf
    }) %>% 
    subset(membership>=.6) %>% 
    (function(mydf) {        
        mydf %>% 
        ggplot(
            mapping=aes(
                x=Var2 %>% factor(
                    levels=coldata_ESCNPCNeu_isoHD2[-c(8,14,21),] %>% subset(grepl('Neuron',Celltype)) %>% arrange(Qlength) %>% row.names,
                    labels=gsub('readsunique.novel.|.bed','',coldata_ESCNPCNeu_isoHD2[-c(8,14,21),] %>% subset(grepl('Neuron',Celltype)) %>% arrange(Qlength,Replicate) %>% row.names)
                ),
                y=value,
                group=Var1,
                color=membership
            )
        ) +
        geom_line(alpha=.8) +
        geom_line(
            data=mydf %>% 
            group_by(Var2,cluster) %>% 
            summarise(meanvalue=value %>% mean),
            mapping=aes(y=meanvalue,group='mean'),
            color='black',size=2
        ) +
    geom_hline(yintercept=0,linetype='dashed',size=.5) +
        scale_color_gradientn(
            colours=colorRampPalette(colors=c('yellow','green','blue','purple','red'))(255),
            limits=c(.6,1),
            na.value=NA
        ) +
        facet_grid(
            facets=cluster~'isoHD: Neuron',
            labeller=labeller(cluster=c('1'='Cluster 1','2'='Cluster 2','3'='Cluster 3','4'='Cluster 4','5'='Cluster 5'))
        ) +
        theme_bw(base_size=16) +
        theme(axis.text.x=element_text(angle=90,hjust=1,vjust=.5),legend.title=element_text(size=12),legend.text=element_text(size=10)) +
        labs(x='',y='Scaled PSI')
    })
) %>% 
ggplot_build %>% 
ggplot_gtable %>% 
with(grobs) %>% 
(function(groblist) {
    groblist[[((groblist %>% sapply(with,name))=='guide-box') %>% which]]
})
`summarise()` regrouping output by 'Var2' (override with `.groups` argument)

In [20]:
options(repr.plot.height=10,repr.plot.width=18)
grid.arrange(
    fuzzyglist[[1]],fuzzyglist[[2]],fuzzyglist[[3]],
    fuzzylegend,
    boxplotglist[[4]],boxplotglist[[5]],
    layout_matrix=matrix(c(1,1,1,2,2,2,3,3,3,4,5,5,5,6,6),nrow=1)
)
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
In [21]:
fuzzygroblist <- fuzzyglist %>% lapply(ggplotGrob)
boxplotgroblist <- boxplotglist %>% lapply(ggplotGrob)

fuzzygroblist <- fuzzygroblist %>% lapply(function(mygrob) {
    mygrob$heights <- boxplotgroblist[[1]]$heights
    mygrob
})
boxplotgroblist <- boxplotgroblist %>% lapply(function(mygrob) {
    mygrob$heights <- boxplotgroblist[[1]]$heights
    mygrob
})
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
Warning message:
“It is deprecated to specify `guide = FALSE` to remove a guide. Please use `guide = "none"` instead.”
In [22]:
#options(repr.plot.height=10,repr.plot.width=18)
grid.arrange(
    fuzzygroblist[[1]],fuzzygroblist[[2]],fuzzygroblist[[3]],
    fuzzylegend,
    boxplotgroblist[[4]],boxplotgroblist[[5]],
    layout_matrix=matrix(c(1,1,1,2,2,2,3,3,3,4,5,5,5,6,6),nrow=1)
)
In [23]:
geneidsuniverse <- hg38samtags %>% 
subset(gene_name %in% (
    filtered.human_ESCNPCNeu_isoHD.leafcutter_highlyused.novel.tsv %>% with(genes) %>% unique %>% strsplit(split=',') %>% unlist %>% unique
)) %>% 
with(gene_id) %>% 
unique

hs_mart <- useEnsembl(biomart='genes',dataset="hsapiens_gene_ensembl",mirror='asia')
genesuniverseBM <- getBM(
    filters = "ensembl_gene_id", 
    attributes = c('ensembl_gene_id','entrezgene_id','hgnc_symbol','description'),
    values = geneidsuniverse,
    mart = hs_mart,
    useCache = FALSE
)
genesuniverseBM$entrezgene_id <- genesuniverseBM$entrezgene_id %>% as.character
entrezgeneidsuniverse <- genesuniverseBM %>% subset(ensembl_gene_id %in% geneidsuniverse) %>% with(entrezgene_id %>% as.character %>% unique)
Batch submitting query [=======>-----------------------]  25% eta:  8s

Batch submitting query [===============>---------------]  50% eta:  6s

Batch submitting query [======================>--------]  75% eta:  3s
                                                                      

3.4. Functional enrichment of Gene Ontology terms

In [24]:
isohdcortexstriatum.mfuzz.list.enrichGOBP <- leafcutterAStable_human_ESCNPCNeu_isoHD.robustintronjunctionsintersectlistlongrenameddf.vastdbevents2.tsv %>% 
subset(JunctionName %in% (
    isohdcortexstriatum.alljunctions.naomit.notstable.membership %>% names %>% 
    subset(isohdcortexstriatum.alljunctions.naomit.notstable.membership>=.6)
)) %>% 
with(ensembl_gene_id) %>% as.character %>% strsplit(split=',') %>% unlist %>% unique %>% sort %>% 
enrichGO(
    OrgDb=org.Hs.eg.db,
    keyType='ENSEMBL',
    ont='BP',
    pvalueCutoff=1,
    qvalueCutoff=1,
    universe=genesuniverseBM$ensembl_gene_id
)
In [25]:
isohdcortexstriatum.mfuzz.list.enrichGOBP2 <- isohdcortexstriatum.mfuzz.list.enrichGOBP %>% 
pairwise_termsim %>% 
simplify(cutoff=0.6,by="p.adjust",select_fun=min)
In [26]:
isohdcortexstriatum.mfuzz.list.enrichGOBP2 %>% 
data.frame %>% 
subset(qvalue<=.1)
A data.frame: 18 × 9
IDDescriptionGeneRatioBgRatiopvaluep.adjustqvaluegeneIDCount
<chr><chr><chr><chr><dbl><dbl><dbl><chr><int>
GO:0043242GO:0043242negative regulation of protein-containing complex disassembly9/252 72/13203 8.694720e-060.013320870.01244551ENSG00000038532/ENSG00000074054/ENSG00000118200/ENSG00000145016/ENSG00000148700/ENSG00000156304/ENSG00000158856/ENSG00000163539/ENSG00000197694 9
GO:0016569GO:0016569covalent chromatin modification 23/252432/132039.021924e-060.013320870.01244551ENSG00000004487/ENSG00000005483/ENSG00000036549/ENSG00000061273/ENSG00000099956/ENSG00000100393/ENSG00000105127/ENSG00000108773/ENSG00000120071/ENSG00000127663/ENSG00000135655/ENSG00000140632/ENSG00000142453/ENSG00000147050/ENSG00000147180/ENSG00000147548/ENSG00000153922/ENSG00000166313/ENSG00000172977/ENSG00000183337/ENSG00000184304/ENSG00000187531/ENSG0000022734523
GO:0010639GO:0010639negative regulation of organelle organization 20/252370/132032.851218e-050.021049110.01966590ENSG00000004487/ENSG00000038532/ENSG00000074054/ENSG00000078018/ENSG00000087470/ENSG00000099956/ENSG00000103197/ENSG00000107960/ENSG00000108773/ENSG00000115317/ENSG00000116670/ENSG00000118200/ENSG00000142459/ENSG00000147601/ENSG00000148700/ENSG00000154473/ENSG00000158856/ENSG00000163539/ENSG00000183337/ENSG00000197694 20
GO:0010975GO:0010975regulation of neuron projection development 23/252485/132035.600666e-050.027564610.02575324ENSG00000004487/ENSG00000078018/ENSG00000087470/ENSG00000091409/ENSG00000100393/ENSG00000101040/ENSG00000107404/ENSG00000118200/ENSG00000127603/ENSG00000142453/ENSG00000142949/ENSG00000147852/ENSG00000148660/ENSG00000160145/ENSG00000163785/ENSG00000166313/ENSG00000169398/ENSG00000174672/ENSG00000184304/ENSG00000197283/ENSG00000197555/ENSG00000198910/ENSG0000021367223
GO:0051893GO:0051893regulation of focal adhesion assembly 7/252 58/13203 1.113041e-040.039021410.03645717ENSG00000065613/ENSG00000074054/ENSG00000115109/ENSG00000127603/ENSG00000158856/ENSG00000163539/ENSG00000169398 7
GO:0090109GO:0090109regulation of cell-substrate junction assembly 7/252 58/13203 1.113041e-040.039021410.03645717ENSG00000065613/ENSG00000074054/ENSG00000115109/ENSG00000127603/ENSG00000158856/ENSG00000163539/ENSG00000169398 7
GO:0031122GO:0031122cytoplasmic microtubule organization 7/252 60/13203 1.384104e-040.039021410.03645717ENSG00000065613/ENSG00000074054/ENSG00000095066/ENSG00000107404/ENSG00000109171/ENSG00000118200/ENSG00000204843 7
GO:0007020GO:0007020microtubule nucleation 5/252 28/13203 1.670246e-040.039021410.03645717ENSG00000074054/ENSG00000109171/ENSG00000141551/ENSG00000163539/ENSG00000204843 5
GO:0007528GO:0007528neuromuscular junction development 6/252 44/13203 1.750543e-040.039021410.03645717ENSG00000107404/ENSG00000151150/ENSG00000160145/ENSG00000167535/ENSG00000198722/ENSG00000204843 6
GO:0010810GO:0010810regulation of cell-substrate adhesion 12/252181/132031.849982e-040.039021410.03645717ENSG00000065613/ENSG00000074054/ENSG00000091409/ENSG00000115109/ENSG00000127603/ENSG00000130402/ENSG00000141503/ENSG00000158856/ENSG00000163539/ENSG00000164620/ENSG00000169398/ENSG00000196924 12
GO:0098732GO:0098732macromolecule deacylation 8/252 103/132037.862923e-040.072560030.06779184ENSG00000061273/ENSG00000073910/ENSG00000100393/ENSG00000105127/ENSG00000143353/ENSG00000184304/ENSG00000187531/ENSG00000196924 8
GO:0016482GO:0016482cytosolic transport 10/252159/132039.513604e-040.078037980.07290982ENSG00000008294/ENSG00000051009/ENSG00000054523/ENSG00000061987/ENSG00000078018/ENSG00000095066/ENSG00000105875/ENSG00000131374/ENSG00000143952/ENSG00000204843 10
GO:0000423GO:0000423mitophagy 4/252 25/13203 1.195359e-030.090510130.08456237ENSG00000038532/ENSG00000103197/ENSG00000115317/ENSG00000161011 4
GO:0051272GO:0051272positive regulation of cellular component movement 19/252461/132031.373791e-030.099003030.09249717ENSG00000008294/ENSG00000037280/ENSG00000054654/ENSG00000060237/ENSG00000061273/ENSG00000074054/ENSG00000078018/ENSG00000087470/ENSG00000091409/ENSG00000105792/ENSG00000115109/ENSG00000130402/ENSG00000131375/ENSG00000158856/ENSG00000161671/ENSG00000163539/ENSG00000169398/ENSG00000184304/ENSG00000196924 19
GO:0051646GO:0051646mitochondrion localization 5/252 44/13203 1.439820e-030.099003030.09249717ENSG00000054523/ENSG00000087470/ENSG00000108773/ENSG00000140983/ENSG00000171109 5
GO:0010632GO:0010632regulation of epithelial cell migration 11/252198/132031.481076e-030.099003030.09249717ENSG00000037280/ENSG00000061273/ENSG00000074054/ENSG00000115109/ENSG00000127603/ENSG00000131375/ENSG00000161671/ENSG00000163539/ENSG00000169398/ENSG00000173482/ENSG00000184304 11
GO:0032259GO:0032259methylation 15/252327/132031.558201e-030.099003030.09249717ENSG00000004487/ENSG00000005483/ENSG00000099956/ENSG00000125835/ENSG00000138050/ENSG00000141644/ENSG00000142453/ENSG00000145388/ENSG00000147050/ENSG00000147180/ENSG00000147548/ENSG00000165792/ENSG00000183337/ENSG00000185917/ENSG00000187531 15
GO:0043551GO:0043551regulation of phosphatidylinositol 3-kinase activity 5/252 45/13203 1.594949e-030.099003030.09249717ENSG00000068078/ENSG00000105875/ENSG00000145016/ENSG00000169398/ENSG00000184304 5
In [27]:
options(repr.plot.height=6,repr.plot.width=8)
g <- isohdcortexstriatum.mfuzz.list.enrichGOBP2 %>% data.frame %>% 
mutate(
    gene.ratio=((gsub('/.*','',GeneRatio) %>% as.numeric)/(gsub('.*/','',GeneRatio) %>% as.numeric)),
    bg.ratio=((gsub('/.*','',BgRatio) %>% as.numeric)/(gsub('.*/','',BgRatio) %>% as.numeric))
) %>% 
mutate(
    enrichment.score=gene.ratio/bg.ratio
) %>% 
head(n=10) %>% 
ggplot(
    mapping=aes(
        x=gene.ratio,
        y=Description %>% factor(levels=Description %>% rev),
        color=p.adjust,
        size=Count
    )
) +
geom_point() +
scale_color_gradient(name='adj. P-val',low='red',high='blue',limits=c(0,.1),breaks=c(0,.05,.1)) +
theme_bw(base_size=14) +
guides(color=guide_colorbar(reverse=TRUE)) +
coord_cartesian(xlim=c(0,.1)) +
labs(x='Gene Ratio',y='Top 10 GO terms') +
theme(axis.text.x=element_text(angle=270,vjust=.5))

g
In [28]:
sessionInfo()
R version 4.0.2 (2020-06-22)
Platform: x86_64-conda_cos6-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /opt/miniconda3/lib/libopenblasp-r0.3.10.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
 [1] tcltk     parallel  stats4    stats     graphics  grDevices utils    
 [8] datasets  methods   base     

other attached packages:
 [1] enrichplot_1.10.1      Mfuzz_2.50.0           DynDoc_1.68.0         
 [4] widgetTools_1.68.0     e1071_1.7-4            reshape2_1.4.4        
 [7] gridExtra_2.3          ggplotify_0.1.0        dendsort_0.3.4        
[10] dendextend_1.15.1      DT_0.16                pheatmap_1.0.12       
[13] clusterProfiler_3.18.0 org.Hs.eg.db_3.12.0    AnnotationDbi_1.52.0  
[16] IRanges_2.24.1         S4Vectors_0.28.1       Biobase_2.48.0        
[19] BiocGenerics_0.36.0    biomaRt_2.46.0         ggplot2_3.3.5         
[22] dplyr_1.0.2           

loaded via a namespace (and not attached):
 [1] fgsea_1.16.0         colorspace_2.0-1     ellipsis_0.3.2      
 [4] class_7.3-17         IRdisplay_0.7.0      qvalue_2.22.0       
 [7] base64enc_0.1-3      farver_2.1.0         graphlayouts_0.7.1  
[10] ggrepel_0.9.1        bit64_4.0.5          fansi_1.0.2         
[13] scatterpie_0.1.5     xml2_1.3.2           splines_4.0.2       
[16] GOSemSim_2.16.1      polyclip_1.10-0      IRkernel_1.1.1      
[19] jsonlite_1.7.2       GO.db_3.12.1         dbplyr_2.0.0        
[22] ggforce_0.3.2        ggvenn_0.1.9         BiocManager_1.30.16 
[25] compiler_4.0.2       httr_1.4.2           rvcheck_0.1.8       
[28] assertthat_0.2.1     Matrix_1.3-4         fastmap_1.1.0       
[31] tweenr_1.0.1         htmltools_0.5.3      prettyunits_1.1.1   
[34] tools_4.0.2          igraph_1.2.5         gtable_0.3.0        
[37] glue_1.4.2           DO.db_2.9            rappdirs_0.3.1      
[40] fastmatch_1.1-0      Rcpp_1.0.7           vctrs_0.3.8         
[43] ggraph_2.0.4         stringr_1.4.0        lifecycle_1.0.1     
[46] XML_3.99-0.3         DOSE_3.16.0          MASS_7.3-55         
[49] scales_1.1.1         tidygraph_1.2.0      hms_0.5.3           
[52] RColorBrewer_1.1-2   curl_4.3             memoise_1.1.0       
[55] downloader_0.4       yulab.utils_0.0.4    stringi_1.4.6       
[58] RSQLite_2.2.1        BiocParallel_1.24.1  repr_1.1.0          
[61] tkWidgets_1.68.0     rlang_0.4.10         pkgconfig_2.0.3     
[64] evaluate_0.15        lattice_0.20-45      purrr_0.3.4         
[67] labeling_0.4.2       htmlwidgets_1.5.3    cowplot_1.1.1       
[70] shadowtext_0.0.7     bit_4.0.4            tidyselect_1.1.1    
[73] plyr_1.8.6           magrittr_2.0.1       R6_2.5.0            
[76] generics_0.1.2       pbdZMQ_0.3-4         DBI_1.1.0           
[79] pillar_1.7.0         withr_2.3.0          tibble_3.1.6        
[82] crayon_1.4.1         uuid_0.1-4           utf8_1.2.2          
[85] BiocFileCache_1.14.0 viridis_0.6.2        progress_1.2.2      
[88] grid_4.0.2           data.table_1.13.6    blob_1.2.1          
[91] digest_0.6.27        tidyr_1.1.2          gridGraphics_0.5-1  
[94] openssl_1.4.3        munsell_0.5.0        viridisLite_0.4.0   
[97] askpass_1.1         
In [29]: